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Abstract 

The bifurcation of an incompressible neo-Hookean thick hyperelastic plate with a ratio 
of thickness to length rj and subject to pure bending is considered within a plane-strain 
framework. The two incremental equilibrium equations corresponding to a nonlinear pre- 
buckling state of strain are reduced to a fourth-order linear eigenproblem that displays a 
multiple turning point. It is found that for < rj < oo the plate experience an Euler-type 
buckling instability which in the limit r] ^ oo degenerates into a surface instability. Sin- 
gular perturbation methods enable us to capture this transition, while direct numerical 
simulations corroborate the analytical results. 

Keywords: hyperelastic plates, incremental equations, turning points, boundary layers. 
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1 Introduction 

The development of compressive stresses in mechanical structures is well known to be responsi- 
ble for Euler-type buckling instabilities. What is less recognised is that such scenarios are likely 
to occur in a number of cases that, apparently, are of a completing different nature. A typical 
example is the phenomenon of stress concentration in perforated thin elastic plates subjected to 
tension. Usually, the holes act as stress concentrators that can be completely or only partially 
surrounded by compressed regions. If the pulling forces are sufficiently strong an out-of-plane 
bending instability is experienced locally near the sites of the holes. A systematic investigation 
of problems of this nature has recently been initiated by Coman et al. [H [21 El S] • 

A second example where Euler-type buckling is indirectly encountered is provided by the 
pure bending of a thin and short elastic tube. The curved configuration adopted by the tube is 
characterised by compressive axial stresses on the concave side, whereas tension will prevail on 
the convex part. Experience shows that a regular instability pattern consisting of many little 
ripples will develop along the former region, eventually leading to the creation of one or several 
kinks that signal the collapse of the tube. 

The stability problem of pure bending in thin-walled tubular structures has a long history 
and there is a vast mechanical engineering literature dealing with various aspects; some of it is 
aptly summarised in the authoritative account of Kyriakides & Corona [S] . On the mathemati- 
cal side, noteworthy contributions in the present context are the works by Seide & Weingarten 
[6] and those by Tovstik et al. (briefiy summarised in [7]). The former investigation is based 
on the Donnell-von Karman buckling equations linearised around a variable-coefficient mem- 
brane state of stress. The resulting boundary value problem was analysed numerically with the 
help of the Galerkin method, and it was found that the circumferential shape of the buckled 
cylinder displays a small dimple on the compressed side. Several versions of the same problem 
have been re-considered in [7j from the point of view of asymptotic analysis. Both works just 
now mentioned made the simplifying assumption that the rippling pattern is the same at every 
point along the axis of the cylinder or, in other words, a solution with separable variables was 
a priori postulated. The assumption is sensible for short tubes (which are fairly stiff), but it is 
inadequate for modelling the collapse in the elasto-plastic regime for moderate lengths, which 
turns out to require a very different approach (cf. [5]). 

Our main aim in the present investigation is to re-visit the pure bending of a rubber block 
deforming in plane strain, a problem that has several points in common with the tube bending 
mentioned above. Unfortunately, the literature in this area has focused mainly on describing 
the deformation itself rather than its potential bifurcations. The typical scenario is outlined in 
Figure [T] the undeformed configuration is shown in the left-hand sketch and is characterised 
by the geometric parameters 2L (length), H (height), and 2 A (thickness); the deformed block 
appears on the right in the same Figure. The plane-strain hypothesis simplifies the problem 
considerably, since one needs deal only with cross-sections (shown shaded) perpendicular to the 
vertical axis of the block and situated sufficiently far away from the lower and upper faces. 

Pioneering work on pure bending stability was carried out by Triantafyllidis [S] , who exam- 
ined incremental bifurcation equations for a couple of piecewise power law constitutive models, 
including a hypoelastic one. He pointed out that the underlying instability mechanism is a 
surface instability similar to that encountered in the plane-strain half-space problems discussed 
by Hill & Hutchinson [Hj or Young [10] . Haughton [TT] performed a similar analysis for hypere- 
lastic materials in a three-dimensional context (neo-Hookean, mostly) and allowed for vertical 
compression as well. The instability was found to be of Euler-type but his interpretation of 
some of the results is wrong (as we shall explain in §3). A novel feature in is the interac- 
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Figure 1: Cylindrical bending of a rubber block; left: reference config- 
uration, right: current configuration. Shaded areas indicate two generic 
cross-sections perpendicular to the vertical axis. 



tion between two different modes of instability, one due to pure bending, the other related to 
compression. Dryburgh & Ogden [12] introduced thin coatings on the curved boundaries of the 
bent block and made comparisons with the uncoated case. Their findings show that, relative to 
the latter case, bifurcation is generally promoted by the presence of surface coating, on either 
or both curved boundaries, that is the bifurcation occurs at smaller strains. The relative sizes 
of the shear moduli for the coating and, respectively the bulk material was found to play an 
important role in describing this phenomenon. 

The finite elasticity works reviewed above share a common feature in that they all deal with 
incompressible materials. So far little is known about the role played by compressibility on the 
bifurcation behaviour in pure bending. The reason might be rooted in the absence from the 
literature of a manageable closed-form expression for the pre-bifurcation deformation. Aron & 
Wang [ini [Hj touched upon issues like existence and uniqueness for bending deformations in 
unconstrained elastic materials, while Timme et al. [15] used Hencky's compressible elasticity 
model to investigate closed-form solutions for cylindrical bending. They succeeded in deriving 
explicit expressions for the bending angle and moment in terms of the circumferential stretches 
on the curved boundaries. The solution is quite involved and it seems unlikely to be useful for 
anything but numerical calculations. Moreover, the particular Hencky elasticity framework is 
restricted by moderate deformations only. 

A critique of bifurcation phenomena in pure bending was given by Gent & Cho [12] who 
pointed out that their experiments did not agree with the theoretical predictions based on the 
surface- instability concept proposed in ^ . In particular, they found that the instability occurs 
for a smaller degree of compression and the block adopts a configuration with a small number of 
sharp creases on the inner surface. To fully explain this observation would require a nonlinear 
post-buckling analysis because the bifurcation involved is probably of subcritical type. We note 
in passing that Gent & Cho's creases are, to a certain extent, very similar to those encountered 
on the curved surface of severely torsioned stocky rubber cylinders [TT]. These phenomena are 
likely to be related to the failure at the boundary of the complementing condition (see [H] and 
the reference therein), and they fall outside the scope of our study. 

With this background in mind, we shall re-consider in the next sections a particular instance 
of the pure bending problems taken up in [TT], [T2]. The aim is to elucidate the nature of the 
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instability and to analyse the mathematical structure of the governing boundary value problem 
when rj = A/L S> 1. To avoid "missing the forest for the trees", the model investigated will 
be confined to incompressible neo-Hookean materials. In §2 these assumptions are showed to 
yield an eigenproblem for a fourth-order partial differential equation with variable coefficients, 
subsequently simplified by seeking a solution with separable variables. Direct numerical sim- 
ulations reveal an Euler-type buckling phenomenon for < rj < oo, but in the limit ?7 — > oo 
this degenerates into a surface instability. Some erroneous interpretations proposed by previous 
investigators (e.g., [H] or [TT]) are also corrected here for the first time. As demonstrated in §4, 
the transition regime between the two different forms of instabilities can be efficiently captured 
by singular perturbation methods. The two contrasting asymptotic methods employed are dis- 
cussed separately in §4.11 (WKB) and, respectively, §4.21 (boundary layers). The former would 
seem to be the most appropriate because the differential equation in question has variable co- 
efficients. However, it transpires that a conventional boundary-layer analysis sheds more light 
and helps us to steer clear from the confusion created by the presence of a multiple turning 
point. Unlike in the recent studies [H [21 El S], here turning points play no role whatsoever 
(the same is true for the related works [25l [27]). The paper concludes with a discussion of the 
results obtained, together with suggestions for further study. 

2 Overview of the model 

Finite pure bending of incompressible hyperelastic materials is discussed in a number of books 
like, for example, [I9] or [20]. To make the paper reasonably self-contained, we summarise some 
of those ideas below. 

The reference configuration of the initially undeformed rectangular cross-section of the hy- 
perelastic block is the region (see Figure [2]) 

■= {(Xi,X2) G I - A < Xi < A, -L < X2 < L} . 

Supposing that the block is bent (symmetrically with respect to the Xi-axis) into a sector 
of circular cylindrical tube, the current configuration of the deformed cross-section is easily 
represented in polar co-ordinates by the domain 

Be := {{r,e) G M X (0,27r] I - n < r < r2, -ooq < < ujq} . 

Rivlin [21] showed that for an incompressible elastic material this type of deformation may be 
described by the mapping 

r={d + 2Xi/u;Y/^, 6 = 00X2, (1) 

where (i is a quantity determined by the particular constitutive law adopted and uj can serve as 
a control parameter as it is related to the angle of bending, ujq := ujL. Since the plate cannot 
be bent into itself, we should require that 

< CJo < TT , (2) 

an assumption used tacitly henceforth. Although the deformation recorded in ([T|) seems to have 
an iconic status among workers in elasticity, it is clear that the kinematics afforded by that 
expression are somewhat restricted. The lines Xi = const, become arcs of the circle r = const., 
while the lines Xi = const, are transformed in lines 6 = const.; in other words, "cross-sections" 
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Figure 2: The undeformed (left) and deformed (right) cross-sections of the 
rubber block shown in Figure [TJ Bending is symmetric with respect to the 
xi-axis so that the two angles marked are congruent and equal to loq; see 
the text for more details. 



perpendicular to the vertical symmetry axis of Br, remain orthogonal to the deformed axis of 
the current configuration Be- This is somewhat at odds with the commonly accepted point of 
view in structural mechanics, according to which pure bending of thick sandwich panels (e.g., 
|24j ) is based on models that allow cross-sections to slide relative to the normal to the deformed 
axis. Nonetheless, the nonlinear mapping ([T]) is still a sensible choice for the type of questions 
we want to answer, at least in a first approximation. 

The bifurcation analysis carried out in this work is based upon linearising the plane-strain 
equations of finite elasticity around the nonlinear pre-buckling deformation ([T]). This approach 
to linearised or incremental bifurcations is well established (see |221 [20], for example), so we 
shall not rehearse it here. Instead, we limit ourselves to pointing out the key steps that lead to 
our eigenproblem. 

We shall assume that the constitutive behaviour of the material is characterised by a strain- 
energy function W = W{\r, Xe), where the principal stretches A,, and Xq are associated with the 
Eulerian principal directions and, respectively, eg. Due to the incompressibility constraint 
these can be written as 

Xr = X~^ and Xg = X := ur , 

which defines the notation A. According to ^20j, the two-dimensional version of the incremental 
equations of equilibrium for incompressible elasticity read 

divs = 0, divu = , (3) 

where u = {u{r, 6),v{r, 9)) is the incremental displacement field, and s denotes the incremental 
nominal stress tensor with components 

°Sij = Lijkihi + pFij - p5ij , i, j ^ {r, 9}. 



Here p is the increment in the Lagrange multiplier p = p{r,9) (the "hydrostatic pressure"). 
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while Fij represent the components of the incremental deformation gradient, 
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Finally, L^j^^i are the components of the fourth-order tensor of instantaneous incremental moduli 
which, in Eulerian principal axes has 15 independent non-zero such components (cf. 



Lijij 



XiXjWij , 



xUxw, - X,W,) 

II % ^ 2^ Aj , 

\ — Aj 

^{Liiii - Liijj + XiWi) if i ^ j, X., 



Xj , 



X^W,, 



with Wi = dW/dXi, Wij = d'^W/dXiXj, and the summation convention does not apply. 
Direct calculations show that the system of equations can be reduced to 



-^1122 + P^r ) + -^1111 + -^2222 ~ 2L1122] U,r 

+ r^{LlUl - Lii22)u,rr +-^2121 (m,^^ -^,9 ) + '"-^2112^'- 



(4a) 
(4b) 
(4c) 



re 



(5) 



rp,g = {r 1^212 + Li212){rV,r +U,g -V) + r'^ L1212V, rr 

+ r{L2ll2 + -^^1122 - L2222)u,r9 ■ (6) 

To avoid overdoing the notation we have used the correspondence r ^ 1 and 6^2 for 
the incremental moduli, and have indicated their derivatives with respect to r by dashes. A 
further simplification is afforded by the incompressibility condition which allows us to deduce 
the existence of a potential = 0(r, 6) such that 

u = -— v = -— (7) 

r 06' dr 

The upshot of this observation is that the two equations can now be combined into a 

single partial differential equation for the potential function. After some routine (but lengthy) 
manipulations, we end up with 

4 

Y.CA<P]=0, (8) 

i=i 

with Cj partial differential operator of the j-th order defined according to 
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The form of the bifurcation equation is vahd for any choice of incompressible hyperelastic 
material but, as it stands, the model is not easily amenable to analytical work. Before further 
simplifications are implemented, we must address the issue of boundary conditions for ([8]). 
The two curved boundaries of Be are taken to be traction-free, a constraint which demands 



ar^(f),rrr -(2/3 + a){(f),ge -r4),ree ) = , for (r, 9) e {n, rs} x (0, 27i] , (9a) 

(j),ee +r^,r -r^^,rr = , for (r, 9) E {n, rs} x (0, 2tt] . (9b) 

These conditions can be obtained as a particular case of the calculations of Dryburgh & Ogden 
|12] . to which the reader is referred for more information. 

Next, we look for separable solutions of the bifurcation equation ([8]) in the form 

(j){r,9) = <^{r)cos{m9) , (10) 

where m G N is the azimuthal mode number related to the number of ripples on the compressed 
side of the rubber block, while $(r) is the infinitesimal amplitude of this cosine rippling pattern. 
Several types of boundary conditions are possible for the straight boundaries of Be- Following 
pTl [12] we consider zero incremental displacement in the radial direction and vanishing normal 
traction. It can be shown that these conditions are satisfied as long as 

nil 

m=— , 11 
ujL 

for some positive n G Z. With this information in hand, all that remains to be done is to choose 
a constitutive model and carry out the simplification of (IHl) with the help of the assumed form 
solution recorded in (|TOl) . 

The bulk material is modelled by a simple neo-Hookean strain-energy function specialised 
to plane-strain elasticity, 

T being the ground-state shear modulus of the material. As shown by Rivlin [21] and subse- 
quently discussed by others [SI [HI [12], for this particular choice of constitutive law the constant 
d in ([T]) is determined by 

^=^(l + V-oV^ 

UJq 

On making use of ( ITOi) in ([8]) results in an ordinary differential equation which, when expressed 
in the non-dimensional variables 

r A 

p:=-, /x:=n7r, ^ ■= J ^ 

can be cast as 

$"" + P(p)<l>'" + Q(p)<l>" + 7^(p)<l>' + 5(p)<l> = 0, on pi<p<p2. (12) 
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Above, the dashes denote derivatives with respect to p and 
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Note that the inner and, respectively, the outer curved surfaces of the current configuration 
become 

Pi,2 = -[(l + V^o')^/^±2r/.;o]'^', (13) 
while the principal stretch in the e^-direction assumes the simple form 

A = cuop. (14) 

The solution of (fT2|) is found subject to the non-dimensional boundary conditions obtained 
from ® via (dUD, 



2 / 2 2 I 

P ( WoP + 
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$ = , for p = pi,2 



$ = , for p = pi,2 



(15a) 
(15b) 



The normal-mode approach has reduced the bifurcation analysis to the study of a standard 
ordinary eigenproblem for $(p) and uoq G (0, vr). While for structural mechanics problems (e.g., 
[23]) this route is free of pitfalls, in finite elasticity it is only deceptively so. The danger is that 
the bent block might develop shear bands or other material instabilities before the compressed 
inner surface starts to wrinkle. Such occurrences are heralded by a loss of ellipticity in the 
partial differential equation ([H]); unfortunately, they remain undetected by (fT^ . Conveniently, 
the use of a neo-Hookean constitutive law precludes any form of material instabilities (note 
that £4 is strongly elliptic in this case). Such exotic effects, however, were accounted for in [8], 
but it was found that the surface instability was always the first to occur. 



3 Numerical experiments 

The stability of the bent rubber block is now investigated numerically, the starting point being 
the eigenproblem fll2|15l) formulate in §2. Our first objective is to find out the dependence of 
the critical bending angle ujq in terms of the aspect ratio t] = A/L. It is expected that an Euler- 
type buckling instability is experienced for < 77 < 00, but in the limit rj ^ 00 this behaviour 
degenerates into a surface instability. Such behaviour is consistent with the results of Dryburgh 
& Ogden [12] and Haughton [11], although it appears that earlier investigators [8] reported only 
surface instabilities for some other choices of constitutive behaviour. Unfortunately, a direct 
comparison with those results is not possible here but, intuitively, one would expect that the 
finite thickness of the block should set a length-scale for the instability pattern. 

A rather peculiar feature of our eigenproblem is the dependence of the mode number m 
in (JTTl) on the bending angle. In order to identify the former quantity for a given rj, we shall 
plot the principal stretch on the curved inner boundary p = pi in terms of this number, the 
critical value of m being that associated with the largest Ai = X{pi,rj;n) when n e N. This 
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Figure 3: Plot of the critical stretch Ai = X{pi) against the mode number m = n/uJo, 
as obtained by direct numerical integration of the eigenproblem (jl2ll5p for a sample 
of aspect ratios rj. The maximum principal stretch for each individual case considered 
is marked by a small circle and corresponds to the points: Pi {rj = 1.0), P2 [r] = 3.0), 
P3 [-q = 5.0), P4 (r/ = 10.0), P5 ij] = 15.0), and P& (r/ = 20.0). These maxima are 
attained for n = 1, except for Pi which corresponds to n ~ 2.43. 



procedure is carried out in Figure [3] where we consider a sample of values for 77 (see the caption 
for details): the horizontal axis records the mode number n/ujQ while the vertical axis shows 
Ai. Strictly speaking, n G N but we shall take this parameter to be a positive real number and 
notice that the eigenvalue of the problem (I12fl5l) . cuq, will depend on this quantity as well as 
on ?7, i.e., ujq = uoQ{ri,n). For rj = 1 we find the curve shown with a continuous line and which 
consists of two sloping parts separated by a peak. Pi (corresponding to n ^ 2.43). Henceforth, 
we shall refer to this curve as Ci. Note that the right-hand part is monotonic decreasing and 
unbounded but here only a segment of that curve is shown. The neutral stability curves for 
the other values oi rj = rjj > 1 considered are Qj := {(Ai(pi; r/j, ra), n/uo) \ n G M+j and they 
all turn out to be part of Ci. The remark made above regarding Ci applies for these curves as 
well. For the sake of clarity in Figure [3] the curves are shifted and shown separately as dashed 
lines, but their top endpoints (P2 ^ Pq) are marked on Ci as well. All such points correspond 
to the choice n = 1. 

The feature illustrated in Figure [3] is generic and not restricted to the particular values of 
7] > 1 chosen. A first observation is that the number of ripples on the compressed side of the 
block increases with the non-dimensional thickness rj. When this latter quantity is reasonably 
large {rj >k, 3) the critical mode number given by (JTTi) always corresponds io n = 1. Thus, 
the behaviour of a very thick block can be understood in two different ways: (z) assuming 
that n = 1 and rj ^ 1 or, conversely, (u) fixing rj = 0(1) and letting n ^ 1. In the former 
case the critical mode number will simply be ir/ujirj, 1), whereas in the latter one the following 
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observation helps: if 771, 772 > are two given, sufficiently large aspect ratios with rji < 772, then 

1 _ m 

Comparing this with (fTTi) assertion (ii) should now be obvious. 

It is instructive to gain some insight into the behaviour of the critical eigenf unctions as- 
sociated with the Pj 's marked on Figure El This information is included in Figure H] where, 
for the sake of brevity, we show the radial and azimuthal displacements only for Pi -h P4, as 
obtained from the two equations in ([71). The localisation of the deformation near the curved 
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Figure 4: The eigenfunctions associated with the critical points Pj in Figure [3l 
(a) Pi (r/ = 1.0), (6) P2 (r/ = 3.0), (c) P3 (j] = 5.0), and (d) P4 (?? = 10.0). In each 
plot the continuous line denotes p~^^{p) {radial displacement), while the dashed 
line is used for ^'{p) {azimuthal displacement). The range for these functions is 
Pi < p < P2 and they are suitably normalised so that their maximum amplitude is 
unity. 



inner surface of the bent block when rj increases is clearly obvious. The stress concentration 
phenomenon revealed by these plots is to be expected because the thicker the rubber block, the 
more difficult is to bend it, that is, the instability will be likely to occur for small values of the 
bending angle. Hence, curvature effects will only be "felt" in the immediate proximity of the 
bent inner surface. In the remaining of the paper we show that this behaviour is ideally suited 
for a singular perturbation analysis. 

At this juncture some remarks on the method used to identify the critical mode number 
are appropriate. At ffist sight, our work in Figure might appear a little awkward. The 
coincidence of the curves Qj {j = 2, ... ,6) with Ci could have been inferred by taking into 
account that the principal stretch Ai is independent of n and depends only on the product uoi] 
- see ( |T3l) and dHl). However, we believe that the longer route taken here has the advantage of 
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clarifying some of the vague statements made by Haughton in [TT] . He misinterpreted the role 
played by n G N in formula ( JTTl) and, in Figure 5 of his paper, he varied both n and r] ending 
up with a wrong statement as to the behaviour of the neutral stability curves for the bending 
problem. For convenience we reproduce that scenario in our Figure Although a different 
formulation of the eigenproblem was used in [TT] (without recourse to any potential function), 
the results we show are the same (as they should be since the height of the rubber block in 
[llj was fairly large, H/A = 10). The only exception is the unusual feature seen in that paper 
for n = 1, which we did not find with our model. Triantafyllidis [8] seems to have committed 
a different error by excluding uq from what he refers to as "wave-number" (see Figure 7 in his 
work). That might explain why he did not find an Euler-type buckling instability. 




Figure 5: A plot of critical values of Ai = X{pi) against undeformed length L/A for 
mode numbers n = 1 10 (see also Figure 5 in reference [11]). The arrow indicates 
the direction of increasing n. 

The response curves shown in Figure are reminiscent of similar situations encountered 
in the buckling of thin- walled structures (e.g., [23]). In that particular type of situation n 
represents the number of half-waves of the instability pattern and plots like the one shown 
above can be used to infer the wavelength of the buckling pattern from knowledge of some 
aspect ratio (related recent work on thin- film instabilities can be found in [H [21 El H] ) • However, 
in the present context such extrapolations appear to provide misleading information for obvious 
reasons. Also, in the limit A/L — > oo the critical principal stretch would have to be equal to 
the value 0.544 predicted by Biot's analysis for a neo-Hookean half-plane in compression (cf. 
[TBI l22] ). This is clearly not the case in Figure [5] but, on the other hand, the earlier results of 
Figure [3] do anticipate this expectation. 
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4 Stress concentration for n ^ 1 

The mathematical structure of the boundary- value problem derived in §2 is akin to a number of 
situations investigated recently in the literature by Fu et al. [2SII2S] and Haughton & Chen |27j . 
Broadly speaking, these authors encountered a particular occurrence of turning points (e.g., 
[28]) or repeated roots in the characteristic equations associated with bifurcation analyses for 
everted cylindrical/spherical shells. It was stated that such special points could aid in detecting 
sites of high-stress concentration within elastic solids. Furthermore, in light of the recent work 
on edge-buckling of thin films [H [21 El S] , it would appear very reasonable to conclude that Fu's 
observation might go a long way towards explaining the localised behaviour seen in Figure |H 
That this is not true we are going to see in §4.21 but before we pursue those issues it is important 
to gain an understanding of the relevance of WKB techniques in the present context. According 
to the previous interpretation of the parameters t] = A/L and n (defined in equation (fTTj) ). the 
bifurcation of a thick rubber block (77 ^ 1) can be understood by taking rj = 1 and allowing 
n ^ 1. This is precisely what we do in the remainder of the paper. 

4.1 WKB approach 

The WKB method is a simple and efficient tool for dealing with variable-coefficient linear 
differential equations containing certain small or large parameters. We shall exploit the presence 
of /i = vm ^ 1 in our eigenproblem to describe the dependence of c^o (or Ai) on this large 
parameter. 

A WKB solution of (fT2|) is sought in the form 

$(p) = Y{p) exp (^/i j^' 5(0 d^^ , (16a) 

Yip) = Fo(p) + -n(p) + ^Y,ip) + ... (16b) 
II p^ 

where 5" = S{p) is one of the roots of the characteristic equation 

V ^oP / 

and Yj{p) (j = 0, 1, ... ) are functions that are to be determined sequentially by substituting 
the ansatz (|T6l) into (JT2l) . and then solving the differential equations obtained by setting to 
zero the coefficients of like powers of p. The above bi-quadratic has four real roots that will be 
labelled 

UqP 

and they lead to a set of linearly independent (approximate) solutions for (|T2l) . Given our 
experience with the direct numerical simulations of ^ it is expected that only the exponentials 
corresponding to S[~2{p) need to be used in order to capture the through-thickness localised 
behaviour. The superscripts "1" and "2" will be used to identify quantities associated with 
these characteristic exponents in ( fT6l) . 

The determinantal equation that follows by imposing the boundary conditions ( fT5|) at p = pi 
on the WKB solutions ^^^^ and has the form 



f/l(pi)^2(pi)-f/2(pi)l^l(pi)=0. 



(17) 
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where 



\/ (p) := $(^)" - -^^^y + , J = 1, 2 . (19) 

When calculating ^^^^ (j = 1,2) we shall ignore terms of order 0(/i~^) and higher in the 
ansatz (116bp . These solutions are fixed by routinely solving a series of non-homogeneous linear 
differential equations. The various multiplicative and additive constants in the expressions of 
those functions can be chosen (without loss of generality) to be unity or, respectively, equal to 
zero. The final results are 



(l-cUoV)^/^' ' 2(l-a;oV 



and thus, 



$(^)(p)^|rj^)(p) + lr/^)(p)|exp(^/i^'5]-)(0d^) , j = l,2. (20) 



It must be noted that Yj''^^'' {j = 0, 1) blow up when p = p = uj^^ G (pi, P2); in the language 
of differential equations this represents a (multiple) turning point of the differential equation 
(|T2l) . Such points are usually defined as those values of the independent variable for which 
some of the roots of the characteristic equation coalesce. For this particular example, both 
'S'i^'*(p) = 'S'2^''(p) and S[~\p) = 5*2"'' (p), i.e., two pairs of roots merge. Strictly speaking, 
the validity of the above formulae for Yj^~^^^ {j = 0, 1) requires |p — p| ^ p~^/^. Although p 
depends on the unknown eigenvalue, our numerical experiments suggest that the turning point 
remains confined to the central part of the interval (pi, P2). 

On making use of (l20l) into (fTTl) . we are able to expand the determinantal equation in 
decreasing integral powers of p 1, 

ro(Ai) + T^{uo, Ai)- + r2(o;o, Ai)^ + ■ ■ ■ = , (21) 
p p 



with 



and 



To{z) := 8z\z^ + lf{z^ - z^ + z + l){z^ + z^ + z - 1) 
Ti{uJo, z) := 2tuo{4:Z^^ + I5z^° + 23z^ + 82^ - 7z'^ - 3) , 



r.(-o, z) := ^^2 + i)2(,._i)4 (6^" - 18-" - 37." " 477^- - lllSz^ 



930z^^ - 600^^° - Az^ - 120z^ - UOz^ + Abz^ + 33) . 
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The solution of (I2T1) yields approximations for both the critical bending angle uq, and the 
principal stretch Ai. For the sake of brevity we record only the final results here 

a;o = Ho + — + ^ + . . . , (22) 

Uq = 0.771844 , Hi = -1.305565 , = 15.39664 , 

and 

Ai = Ao + ^ + ^ + . . . , (23) 

Ao = 0.543689 , Ai = 0.385922 , Aa = -4.184333 . 

As one would expect, Aq ~ 0.544 represents the critical value of the principal stretch for the 
surface instability of a compressed neo-Hookean half-space (cf. [22]); the next-order corrections 
in formula fl^Hl) account for the finite size of the rubber block. To assess the usefulness of the 
two asymptotic results fl2^ and fl2Hl) . a set of comparisons with direct numerical simulations is 
recorded in Table [H The agreement is excellent for both cuq and Ai; in particular, we find that 
the relative accuracy (RA) associated with c^o ranges between 1.4% (n = 7) and 0.8% (n = 20). 
The approximation of Ai is even better, for RA is at most 0.4% (n = 7) in all cases considered. 

The WKB analysis laid out above has the advantage of producing a robust approximation 
for Uq (or Ai) with minimum effort. The presence of the turning point, however, is worrying 
because it tends to obscure the true nature of the localised behaviour exhibited by (fT2ll . It is not 
immediately clear whether such behaviour has anything to do with the turning point and, thus, 
a change of tack is imperative. It will shortly become obvious that conventional boundary- 
layer techniques are better suited for understanding the underlying mathematical structure 
responsible for the scenario depicted in Figure HI The details of that particular approach are 
highlighted next. 

4.2 Boundary-layer analysis 

To begin, we introduce the stretched variable X = 0(1) such that p = pi + Xp~^, and look for 
solutions of (1121) with 

W{X) = WoiX) + WiiX)- + W2{X)\ + . . . , (24a) 

/i /i^ 

o'o = ^^0 + — + ^ + • • • , 24b 

pi = Ao + ^ + ^ + . . . . (24c) 

Although ( 124b[) and f l24cl) are not independent, it helps to expand pi in the form suggested 
here. Of course, when solving the governing equations for the coefficients Wj{X) (j = 0, 1, . . . ) 
one has to remember formula (fT3l) and replace the A^'s with their expressions in terms of the 
^- (j = 0,l,...,/c). 

On substituting into (fT^ we find a hierarchy of differential equations 

£blto=i:i:4"^ (25) 

i=0 j=l 
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Table 1: Comparisons between direct numerical simulations of the eigenproblem 
(|12|15|) and the asymptotic results recorded in the formulae and ([^5]) . The 
bending angle and the azimuthal stretch on the inner boundary associated with the 
former set of values are denoted by and, respectively, Aj'"™. The corresponding 
asymptotic quantities are identified as cJq^-^ and A^^'^. 



n 


UJ 


num 





asy 



A 


num 
1 


A 


asy 
1 


7 


0. 


,744313103 


0. 


,733652725 


0. 


,555301084 


0, 


,552585690 


8 


0. 


,744272309 


0. 


,736778030 


0. 


,554325988 


0, 


,552419947 


9 


0. 


,744928414 


0. 


,739453964 


0. 


,553494699 


0, 


,552104104 


10 


0. 


,745886633 


0. 


,741762924 


0. 


,552780068 


0, 


,551733662 


11 


0. 


,746957132 


0. 


,743772049 


0. 


,552160229 


0, 


,551352710 


12 


0. 


,748046186 


0. 


,745533828 


0. 


,551618218 


0, 


,550981720 


13 


0. 


,749107553 


0. 


,747090530 


0. 


,551140475 


0, 


,550629796 


14 


0. 


,750119338 


0. 


,748474960 


0. 


,550716526 


0, 


,550300415 


15 


0. 


,751072409 


0, 


,749714126 


0. 


,550337796 


0, 


,549994245 


16 


0. 


,751964384 


0. 


,750829238 


0. 


,549997574 


0, 


,549710574 


17 


0. 


,752796399 


0. 


,751837997 


0, 


,549690282 


0, 


,549448050 


18 


0. 


,753571374 


0. 


,752754745 


0. 


,549411416 


0, 


,549205075 


19 


0. 


,754293016 


0. 


,753591537 


0. 


,549157200 


0, 


,548980000 


20 


0. 


,754965302 


0. 


,754358139 


0. 


,548924583 


0, 


,548771235 



in which 

/ ^2 1 \ d 

'0 + 72 



is the boundary-layer (BL) differential operator and Co '■= '^o^o- The quantities A[j^ = aIj\x) 

will be defined as we go along, with the convention that A^-^"* =0. In contrast to the WKB 
analysis the general solution of each one of the equations in ( !25|) is trivially found, for £bl has 
constant coefficients. 

The equations (1251) are solved subject to two types of boundary conditions. The first set is 
obtained from (|T5l) with the help of the ansatz (1241) . and can be cast in the general form 

k—l 

^1 m = E E 4^ ^ ' for X = , (26a) 

i=0 j=0,l 

k—l 

^2[l^.] = EE^fT§' forX = 0, (26b) 



where 



i=0 j,=0,l 

^ d^ f^^ 2\ d . ^ d^ 1 



the remark made for the A^j^ 's applies to the boundary coefficients B^j^ and C^j^ as well. 

The second set of boundary conditions is motivated by the numerical experiments illustrated 
in Figure H] and involves the requirement that 

^^0 as X-. 00, (27) 
dX^ ^ ' 
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for z > 0, j = 0, . . . , 3. It might also be worth pointing out that the solution of equation 
(|T2l) is exponentially small in the outer layer (cf. §4.ip . so that (!27j) can be viewed as matching 
conditions between the inner and the outer solutions. 

The leading order problem for Wq{X) is homogeneous and consists of the differential equa- 
tion fl2^ for k = 0, together with the boundary conditions fl2U|) . Rejecting the exponentially 
growing contributions and imposing the normalisation condition VFo(-^ = 0) = 1, it follows 
that ^ 

Wo{X) = exp(-CoX) - exp(-X/Co) • (28) 

When this function is substituted into the boundary conditions, one obtains an algebraic equa- 
tion, 

Co' + 2Co'-4Co' + l = 0, 

whose only acceptable solution is Co ~ 0.543689. The result is identical to Aq found in §4.11 and 
the same turns out to be true for Qq in f l24b[) . 

The next order problem corresponds to taking /c = 1 in ( 125|) and f l26l) . The coefficients that 
appear in these equations are 



.(1) 

^01 




.(1) 

^02 


= -^(Aofii - 
So 


hAl^]o + ^^oX)(Co'-l), Ag^:=^ 

^0 


r(1) 
^00 


:=-^(Co' + 2). 

So 


^01 


:= ^{AoQi - 
So 


|-Aif]o)(Co-2), 


^(1) 


:= -|(Aol^i + Air]o), 
So 


•-^oi 


- 2^03 • 





The first-order correction Qi in the expansion (]24bl) of the eigenvalue Uq is recovered by enforc- 
ing the Fredholm solvability condition on the non-homogeneous problem satisfied by Wi{X). 
The task is simplified by the observation that the homogeneous problem for Wo{X) is self- 
adjoint. Standard calculations show that this constraint amounts to 

C {^(0)}' + {cg> - Bg>} wm'^m - {wmr 

= (29) 

Notice that the integral on the right-hand side in 0291) is evaluated analytically, and thus the 
solvability condition will reduce to a linear equation in Qi. The solution Qi ^ —1.305562 is, 
for all practical purposes, identical to the WKB result obtained earlier. 

The pattern of the boundary- layer approach for the pure bending problem is now clear: at 
each step one will have to impose a solvability condition for finding Qj that features in fl24bp . 
and then solve (exactly) a non-homogeneous fourth-order boundary-value problem in order to 
get Wj{X). The algebraic manipulations become increasingly unwieldy as we move to further 
orders, but symbolic algebra packages help considerably. We have imposed the solvability 
condition for the W^2-problem and found that the value of Q2 predicted agrees with Q2 to 
within five significant digits; for completeness, the coefficients needed to set up that problem 
are recorded below 

.(2) ._ ^ (M _ o\ a{2) ._ .(1) .(2) ._ .(1) 

^11 ^3 l,So ' ^12 ^02 5 ^13 ^03 ' 

So 
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4? := 74 [^o(^ + Ai)(Co' + 9) + 2(^1 (Co' + 3)Co] , 4? := "^2 + 
so ^0 



4J := ^\ 2Aono(Aon2 + A2(]o)(Co - 1) + [^o(^ + ^i)' + A^fi^] (^4 ^ 3) 



+ 4ni(X + AO(C + l)Co-3fi^Co' , 



4? := -4 [Ain^(^4 _ ^ 2fii(Co' - 2)Co] 



:= [2(Aofi2 + A2l]o)(Co - 2)Co + (Agfi? + Alnl){C^ + 6) + 4Ail]i(Co' + 2)Co] 



r{2) ._ o(l) o(2) ._ o(l) 

C!S = -4 [3(A2fi2 + a2^2) _ 2(,{Q,A2 + Aon2 - ^A.Q^j^ , 
So 

(2) _ _ Ai (2) _ (1) (2) _ ^(1) 

01 ~ A 2 ' ^10 ~ '-'00 ' '-"11 ~ '-'01 • 
^0 

It might appear that our boundary- layer approach is a by-product of adopting the simple neo- 
Hookean form for the constitutive response of the bulk material. However, this impression is 
only apparent, for choosing a strain-energy function of the form 

the same mathematical structure persists. In this case the interpretation of the asymptotic 
results become more difficult because now loss of elliptiticy will be encountered for some value 
of g 7^ 2; such issues will be discussed elsewhere [29]. 



5 Concluding remarks 

We have re-examined the bifurcations in cylindrical bending of a thick rubber block under 
the assumption of plane-strain deformation. The constitutive behaviour was taken to be that 
of a neo-Hookean incompressible solid, the reason for this being twofold: we wanted to (a) 
exclude any material instabilities and, (6) simplify our equations as much as possible. The 
outcome turned out to be a simple fourth-order eigenproblem with variable coefficients. Direct 
numerical simulations and singular perturbation methods were employed to unravel the origins 
of the rippling pattern triggered on the compressed face of the block, when the bending angle 
is sufficiently large. It has been shown in ^ that previous investigators [HI [H] misinterpreted 
the definition of the so-called "mode number" and thus made several erroneous statements. In 
particular, we want to re-iterate here that blocks of a large but finite thickness will always 
experience an Euler-type buckling instability with a well defined number of ripples. It is only 
in the limit of an infinitely thick block that one finds the degenerate surface instability. Our 
results deal with the neo-Hookean material, but it is believed that the above statement remains 
valid in other well. It seems quite unlikely that the choice of constitutive law would 

have any dramatic repercussions on the overall arguments presented in this paper. 

Our work has also demonstrated that the use of WKB methods in the context of incre- 
mental elasticity is unnecessary and even misleading. The multiple turning point featuring 
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in the present eigenproblem has httle to do with the tendency of the ripphng deformation to 
confine itself near the inner curved surface of the block. Interestingly though, a simple-minded 
boundary-layer analysis was able to expose the nature of the localisation in a trivial way. In 
spite of the original complexity of the problem, we found that if the rubber block is sufficiently 
thick, its possible bifurcations from the cylindrical configuration are governed by constant- 
coefficient differential equations - easily solvable in closed form. This has important overall 
implications since preliminary calculations indicate that the problems taken up in [251 [26t [27]) 
are amenable to a similar boundary- layer analysis. We shall report the corresponding details 
in a forthcoming study [22] ■ 
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